Negative differential resistance in molecular junctions: application to graphene ribbon 
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Using self-consistent calculations based on Non-Equilibrium Green's Function (NEGF) formal- 
ism, the origin of negative differential resistance (NDR) in molecular junctions and quantum wires 
is investigated. Coupling of the molecule to electrodes becomes asymmetric at high bias due to 
asymmetry between its highest occupied molecular orbital (HOMO) and lowest unoccupied molec- 
ular orbital (LUMO) levels. This causes appearance of an asymmetric potential profile due to a 
depletion of charge and reduction of screening near the source electrode. With increasing bias, this 
sharp potential drop leads to an enhanced localization of the HOMO and LUMO states in different 
parts of the system. The reduction in overlap, caused by localization, results in a significant reduc- 
tion in the transmission coefficient and current with increasing bias. An atomic chain connected 
to two Graphene ribbons was investigated to illustrate these effects. For a chain substituting a 
molecule, an even-odd effect is also observed in the NDR characteristics. 



c 

o 
o 



> 
o 

OO 

o\ 
cn 

in 
o 

00 

O 



X 



PACS numbers: 73.23.-b,73.63.-b 



I. INTRODUCTION 



Negative differential resistance (NDR) was first ob- 
served by Esaki in diodes 1 -, where occupied states on one 
side become aligned with the gap of other side as the 
voltage is increased. Current reduction also occurs when 
the position of the resonant states of the molecule move 
within the gap of one of the contacts^ 3 - as in resonant tun- 
neling diodes. In metallic carbon nanotube junctions 4 -, it 
was found that the reduction of the current is due to a 
mismatch in the symmetry of the incoming and outgo- 
ing wavcfunctions of the same energy. Another work^ on 
the I-V characteristic of CoPc on gold has also associated 
the NDR effect with lack of orbital matching between Ni 
tip and Co atom. Another origin was explained in STM 
measurements^ 7 -. In this case, narrow peaks in the local 
density of states (LDOS) of an atomic scale tip sweep 
past the LDOS of an adsorbed molecule as the bias volt- 
age is increased. 

More recently, more instances of NDR were 
observed£i£i2 or predicte d 10 ' 11 ! 12 ' 13 ^ 4 ' 15 in molecu- 
lar devices. In the case of potential barriers in 2D 
Graphene sheets 1 ^., the effect was due to the linear 
dispersion of (massless Dirac) electrons which show a 
gap in their transmission across the barrier. In Ref. [ll| it 
was due to the presence of Van Hove singularities in the 
DOS of the ID electrodes regardless of the type of the 
contact. This latter explanation is related and similar to 
that of Refs.[fj|7||l2j which involves sharp features in the 
LDOS. In these cases, however, the general conditions 
necessary for the observation of the effect were not 
clearly elucidated. Sharp features in the LDOS can lead 
to ND R 12 ' 16 , but it is not a sufficient condition for the 
observation of NDR, as a reduction in spatial overlap of 
those states is also needed. 

The current in nanoscale devices is given by the Lan- 
dauer formula (see eq. I15|) which involves the transmis- 



sion coefficient given by the product of the local den- 
sity of states (LDOS) of the left and right electrodes by 
the off-diagonal matrix elements of the Green's function 
(GF) connecting the left electrode to the right one (see 
eq |16l) . A reduction in the current is caused by a lower- 
ing of either term in the transmission coefficient. While 
NDR in some devices is caused by a lowering of the ma- 
trix clement of the GF 17 , in some other cases it is caused 
by a reduction in the product LDOS within the energy 
integration window-£&£ i 12 ' 14 i 16 . 

In this paper we explain the reason for occurrence of 
sharp features in LDOS, and also emphasize that charg- 
ing effects play an enhancing role in producing NDR in 
the I-V characteristics of nano-j unctions. A large bias 
causes charge depletion, an asymmetric potential profile, 
and asymmetric coupling even in a symmetric structure, 
resulting in a stronger localization of states on different 
parts of the system, thereby reducing transmission and 
current. 

We consider an atomic carbon chain between two 
graphene tips as a nano-junction (Fig. [I), albeit all re- 
sults are generalizable to other types of nano-j unctions. 
Weak contacts between tips and the chain/molecule 
which usually occur in experiments involving break or 
molecular junctions, are necessary for causing localized 
states within the molecular region and observation of 
NDR. So, we adopt a model in which hoppings to leads 
are smaller than intramolecular or intralead hoppings. 
We claim that in molecular junctions where NDR is ob- 
served, localization of electronic states within the bias 
energy window is the dominant cause of reduction in 
current. The weak bond can play the role of a barrier 
to localize states within or near the molecule. The pur- 
pose of our model is not to make quantitative predictions, 
but just to illustrate the NDR mechanism using a simple 
enough model. Given the small size of contact we assume 
that transport at high bias is mostly coherent and dissi- 
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FIG. 1: (Color online) Two Graphene tips connected with a 
weakened bond to a carbon chain. Tips have sharp structures. 
The weakened bond is considered to be 0.3 of the normal 
hopping of the C-C bond. The central interacting region is 
shown with the dashed rectangle. 

pation due to electron-phonon interactions occurs mainly 
in the drain. 

After presenting Hamiltonian of the system in section 

II. we will introduce the formalism and method used to 
handle the electrostatics of the problem in the section 

III. In the Appendix, electrostatic potential calculated 
by this method is compared with two other methods. 
We are going over the general formalism used for the cal- 
culation of non-linear transport characteristics in section 

IV. The responsible for current reduction in an atomic 
chain between two graphene tips which is known to be 
localization of states induced by charging effects will be 
presented in the section V. 

II. MODEL 

The single electron Hamiltonian of the central system 
(C) including the molecule is 

He = J> + uT l + W^cta + J2 Mcj + *ct) (1) 

ieC <ij> 

where cj and Ci are respectively the electron creation 
and annihilation operators on site i of C, and t is the 
hopping energy between nearest neighbor atoms. One 
7r orbital per site is considered for this system. Un- 
der an applied bias, the solution to Poisson's equation 
is the sum of the solution to Laplace with symmetric 
boundary conditions on the electrodes V(z = 0) = — V/2 
and V(z = L) = V/2 (this is denoted by uf xt ), and 
the solution to Poisson with boundary condition V(z = 
0) = V(z = L) = at both ends (this is denoted by 
Wi = J2j Vij$ n j)- The sum u ext + W clearly satisfies 
Poisson equation and the proper boundary conditions. 
Here Vij is the electrostatic Green's function calculated 
by the method of images, and Snj — rij — is the change 
in the self-consistent charge rij from its initial equilibrium 
zero-bias value. 

It should be noted that parts of electrodes (here also 
called as " tips" ) have been incorporated inside the inter- 
acting central region as there is always some potential 



drop beyond the contact of the electrodes with the cen- 
tral "molecule". 



III. ELECTROSTATIC GREEN'S FUNCTION 

The electrostatic potential is determined by both the 
direct interaction of electrons with each other and the in- 
direct one via image charges. The image charges induced 
by electrons within the electrodes, strongly depend on 
the spatial configuration of the electrodes and the con- 
tact atoms. For the simplicity of calculations, it is usual 
to consider the electrodes as two infinite planes perpen- 
dicular to the molecule^. These planes are located on 
the contacts. 

It is supposed that electrodes are a perfect metal with 
good screening properties, and that at their boundary 
the potential can be considered as a constant so that 
Dirichlct boundary conditions can be applied there. In 
this case, the potential drop occurs within the central 
part of the sample, which we call "molecule" , although, 
strictly speaking this central region is taken to be larger 
than the molecule itself as there is always some potential 
drop at the contact of the electrodes with the central 
"molecule" . 

It should be mentioned that the 3-dimensional Poisson 
equation needs to be solved in order to find the correct 
potential profile along the molecule. Indeed the electric 
field lines are not necessarily straight lines, and a ID so- 
lution would be incorrect. So the Coulomb Kernel needs 
to be more like the 3-dimensional l/\r — r'\ rather than 
the 1-dimcnsional \r — r'\. 

As there is a finite charging energy when the two elec- 
trons are on the same site, there should be no diver- 
gence in the kernel, and the onsite Coulomb repulsion 
has been modeled by the so-called "Hubbard" parameter 
Uh, which could also contain exchange and correlation 
effects if appropriately chosen. However, image charges 
potential lowers the potential on one site from its initial 
value Uh- 

In this article, the Ohno-Klopmann (OK) modet^ has 
been adopted for the Coulombic function U : 



U(n, r} ) = = (2) 

y/\ n - f 3 |2 +U H 2 

It has the correct limits for both large and small inter- 
particle distance ri—rj. It has the advantage of including 
onsite correlations through the Hubbard-likc parameter 
U H - 

In the literature^, there exists an exact Dirichlet 
Green's function for a point charge or a distribution of 
charges between parallel conducting planes held at zero 
potential. The planes are located at z=0 and z=L. Using 
this Green function, we present the following exact form 
which is appropriate for the kernel of Ohno-Klopmann 
model (Eq.@). 
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V(x,y,z;x ,y ,z ) = 2 / dkJ (ak)f(k,z < ,z > ) (3) 
Jq 



where 



f(k,z <1 z > ) 



sinh(fcz<) s'mh(k(L — z>)) 
sinh(fcL) 



a = \J (x - x') 2 + {y - y') 2 + U H 2 



(4) 



(5) 



The asymptotic behavior of the function f(k) in Eq.j4|) 
is as follows: 



lim f(k) -> 0.5e- fc(z >~ z < ) 



0.5 




z< = Z> 



(6) 



Moreover, f(k) goes to zero when k — > 0. Since at z < = 
2>, the function of f(k) will be a constant for k 3> l/z<, 
the integration with infinite range can be converted to a 
limited range integration . 



V( Z< 



*>) = -- 



(1 - 2f(k))J (ak)dk (7) 



where /(&o) = 0.5. The value of ko in nanotubes and 
graphenes used here is about 100. This value depends on 
the distances between atoms of a molecule and also on the 
distances between two boundary planes (L). In case of on 
site electrostatic potential (x = x ; y = y ), the first term 
of Eq.Q is the Hubbard energy. However, a subtraction 
term which depends on the distances between atoms and 
L, lowers the Hubbard energy from Uh- This term is 
the image charges potential which was considered in the 
variational method, too (Appendix. B). The value of the 
semi-empirical Hubbard term for carbon^, is about 10 
eV=0.37 a.u. So U^ 1 = 2.72 whereas the typical bond 
length is of the order of 1.4 A=2.6 a.u. 

In the Appendix (A,B), we compare this method 
(namely the exact method) with two other methods so- 
called the variational and image charges method. 



IV. CALCULATION OF CHARGE AND 
CURRENT 



,22,23 



The charge is obtained using the NEGF formalisr 
The electrodes electrochemical potentials and the fermi 
functions are shown by (J-l.r and /l._r, respectively. The 
retarded Green's function matrix is: 



G{Z) = [ZI-H-^ L -^ R ]- 1 



(8) 



where Z = E + ir\ is a complex variable whose real 
part is energy and r\ — » + . "I" is the unit matrix. H is 



the molecule Hamiltonian defined by Eq.{T]) in the tight- 
binding approach. , R are the retarded self-energies 
arising from scattering by the left /right semi-infinite elec- 
trodes. These self-energies depend on space configuration 
of the electrodes and the quality of the electrode-molecule 
couplings. We have to obtain the surface Green's func- 
tion of semi-infinite electrodes g P (E) in order to deter- 
mine the self-energy. The Lopez-Sancho's method^ has 
been used to calculate the surface Green's function. The 
retarded self-energies are given by: 



K - r T v gl(E)r v p ee L/R 



(9) 



where t p is the coupling matrix between the electrodes 
and the molecule^ 2 .. Since the hopping terms are short- 
ranged, most elements of the coupling matrix are zero. 
Broadening of the molecule energy levels due to attach- 
ment to the electrodes is related to the self-energies as: 



T p =i\ET p - E£] = 2tttJLDOS( P , E)t p (10) 

Note that the broadenings are proportional to the local 
density of states at the connecting sites to the electrodes. 
It should be noted that in this paper transport is assumed 
to be coherent. The charge density is the sum of two sep- 
arate parts coming from equilibrium and non-equilibrium 
charges. Since the voltage division is symmetric on the 
electrodes, the equilibrium charge n eq is calculated from 
the retarded Green's function as: 



Mo-V/2 



Im[G r u (E)]dE 



(11) 



where /xq = fiR = fiL- The initial charge n° is cal- 
culated by the above integration in zero bias. In the 
non-equilibrium situation, the lesser Green's function 
—iG < (E) represents the occupation number in the pres- 
ence of the two electrodes subject to a bias. The non- 
equilibrium charge n non _ eq is determined in the presence 
of an external bias V. 



non—eq 



2vr 



Mo+V/2 



Mo-V/2 



[-iG<(E)]dE (12) 



It can be simply shown that in the coherent regime 
the lesser Green's function is determined by the retarded 
Green's function (Eq.®). 



non — eq 



2tt 



a + V/2 



ao-V/2 



[G r (T L f L +T R f R )G a } H dE (13) 



where f p = 1 / [1 +exp( E r ^ )] shows the fermi function 
of the electrodes. Finally, both parts of the charge are 
summed to give the total charge: 
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non — eq 



(14) 



Since the molecular Hamiltonian itself depends on the 
electron density, one needs to do a self-consistent process. 
The self-consistent algorithm follows these steps. At the 
first step, the left and right self-energies in Eq.© are 
calculated once before the self-consistent loop. In the 
second step, the Hamiltonian is set using a guess input 
charge. 

The calculation of charges in Eqs. (|ll|13| is a hard step 
as it needs to be well converged. The new and old charges 
can be mixed with each other by using linear mixing or 
Broyden's method^. Using the mixed charge, this pro- 
cess will start from the first step and continue till con- 
vergence is achieved. Finally, having the self-consistent 
charge and potential profiles, the current passing through 
the molecule is calculated by the Landauer formula—. 



2e r°° 

I(V) = T y dET(E,V)[f R (E)-f L (E)} 



2e 

Ti 

2e 
~h 



— oo 
H0+V/2 

HO-V/2 



dET(E, V) 



(15) 



where the second expression is written for zero temper- 
ature. The transmission coefficient T(E, V) is defined 
as: 

T = Ti-[G r T R G a r L ] oc LDOS(L)LDOS(i?)|G Lfl | 2 (16) 

The integral evaluation for charge density in 
Eqs. pill3p has to reach a reasonable accuracy. The 
speed of the convergence process depends strongly on the 
accuracy of the integration process. For weak couplings, 
the van Hove singularities in the density of states (DOS) 
will make it tremendously difficult to integrate the DOS 
along the real axis with desired accuracy. Indeed, the 
singularities arising from the poles of the Green's func- 
tion are close to the real axis. However, in the complex 
energy plane, the DOS along the complex contour away 
from the real axis is very smooth^. The resultant for- 
mula for a contour integration of the equilibrium charge 
is: 



eg 



P 



1" Jo 



(17) 



Mo - V/2 - E 



mm MO V/2 -f- E m i n , . 

' z a = o ( 18 ) 



where E m i n is chosen to be lower than the lowest eigen- 
value of H. 



V. RESULTS 

Fig-© shows the NDR phenomenon in the I-V curves 
of odd and even length chains located between two 
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FIG. 2: (Color online) I-V curves for two graphene tips con- 
nected to the chain with a) even and b) odd chains. The 
hopping of the weakened bond is 0.3 times that of the intra- 
chain hopping (t). In case of "0 atom" where two tips are 
facing each other with no chain in between, the hopping in- 
tegral is equal to 0.1 t. NDR is also observed in two (5,5) 
capped nanotubes (NTs) with a 6 atoms chain in between. 
Current through NT system is 50 times larger than shown. 



graphene tips. To show that NDR is also obtained with 
gapless leads, we have also made calculations for a (5,5) 
carbon nanotube and still observed a reduction in the 
current due to localization of states at the caps of the 
tubes at high bias. Details will be reported elsewhere. 
The NDR threshold voltage for odd length chains is 
higher than that for even chains. The origin of this dif- 
ference can be traced back to the distance of those levels 
which play a role in the observation of NDR from the 
Fermi level. For odd chains, the state at the Fermi level is 
an extended state over the length of the chair*21, whereas 
even chains have a gap at the Fermi level. Therefore typ- 
ically a twice larger bias is needed to observe NDR in 
odd chains compared to even chains of similar length. 

To understand the origin of NDR in this system, in 
Fig. [3] we compare the transmission coefficients at the 
current peak and valley voltages. As one can see from 
the figure, there is a large reduction in the transmission 
of the resonant states when the bias is increased. We will 
show that the reason for this can be traced back to a loss 
of LDOS overlap of the left contact with the right one. 

In Fig. |H the electrostatic potential energy and trans- 
ferred charge (5n = n — hq) profiles are plotted for differ- 
ent biases. These distributions are obtained for a small 
voltage (0.2 V) and voltages of the peak and valley of the 
current. In the linear regime, potential is nearly sym- 
metric. However, by increasing the bias, some charge 
is depleted from the source, thereby weakening the ef- 
fect of screening and enhancing the potential drop fur- 
ther at the source. The asymmetry in the voltage drop 
can be understood in the following way. The transferred 
charge between electrodes and molecule depends on the 
quantum capacitance of the molecule. Quantum capaci- 
tance increases with the surface density of states at the 
source or drain electrochemical potentials. Fig. [5] (a,b) 
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FIG. 3: (Color online) Transmission coefficient through a 
chain connected to two graphene tips. Transmission is plot- 
ted for a chain with a) 4 atoms (even chain) b) 5 atoms 
(odd chain) at peak and valley voltages. Vertical dashed, 
dotted and dash-dotted lines identify the Fermi level, inte- 
gration windows at current-peak voltage and current-valley 
voltage, respectively. A large reduction in the transmission 
can be noticed at higher voltage. Transmission through other 
odd/even chains have similar features. 
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FIG. 4: (Color online) Potential and transferred charge 
(rii — n°) for a chain with 4 atoms between two Graphene 
tips. Profiles for three voltages are plotted; a small voltage, 
and voltages which correspond to the current peak and val- 
ley. Potential and charge has been averaged on each Graphene 
layer. Source is on the right and drain on the left. 



shows that LDOS(Ep + V/2) on the surface layer of the 
source side is much smaller than LDOS(Ep — V/2) on 
the drain side. Due to its capacitive coupling with the 
drain, one state (see Fig. EJa)) which is localized on the 
drain side of the molecule is pinned at Ep — V/2. So 
LDOS(Ep — V/2) remains large as the bias is increased, 
while LDOS(Ep+V/2) gradually decreases when the res- 
onant states in Fig. [5Jb) move away from Ep + V/2. This 
asymmetry in LDOS translates into an asymmetry in the 
couplings of the central region to leads, even though there 
geometric symmetry is enforced. On the side with weaker 
coupling (source side in our case) screening would be less 
effective and potential drop more pronounced. Therefore 
essentially the asymmetry at large biases develops due 
to the asymmetry in the distribution of molecular states 
around the Fermi level. This phenomenon is expected 
to be universal in molecular double junctions with weak 




FIG. 5: (Color online) Surface density of states of a) the 
source and b) drain electrodes on z~0 and z=L shows an 
asymmetric coupling to the chain. In this case, the chain 
contains 4 atoms between the two Graphene tips, c) Total 
charge depletion (Sn) of the central region versus applied bias. 



couplings. Another consequence of the effective weaken- 
ing of the couplings to the leads is the sharpening of the 
molecular states. States near the weak coupling will have 
narrower peaks at high bias. This is a signature of their 
enhanced localization. 

The strong reduction in the transmission arises from 
the localization phenomenon which occurs due to the 
sharp linear potential drop near the source tip. The 
onsite energies are most negative on the left side while 
they are most positive on the right side of the source tip 
(atoms located on 20 and 25 A on Fig. 2]). Therefore 
the LDOS of the left side atoms is large at low ener- 
gies, whereas that of the right side atoms becomes large 
at high energies. This situation is very similar to an 
ionic bond with a large onsite energy difference. The 
bonding and antibonding eigenstates become farther sep- 
arated (compared to when onsite energies were equal), 
and this causes transfer of charge to the low energy site, 
and enhanced localization of orbitals on the sites due to 
the large electric field present. 

The upper half of each curve in Fig. [5] shows LDOS on 
the left and right side atoms of the source tip at voltages 
of the peak and valley of current. It was checked from the 
LDOS data that states with higher energies become lo- 
calized on site labeled by 25 A, while lower energy states 
become localized on the left atoms of the source tip (site 
labeled by 20 A). Therefore the product LDOS at these 
two sites is reduced with increasing bias, due to a re- 
duced overlap, leading to a decrease in T(E) according 
to eqUU In Fig. [S]and for a bias voltage of 0.55 V, strong 
localization occurs at E = — 4.25eV where the transmis- 
sion is also reduced. The lower half of the curves in Fig. 
[6] shows that the transmission closely follows the prod- 
uct of LDOS of the left and right atoms (atoms located 
on 20 A and 25 A°) of the source tip. By increasing the 
bias from the current-peak to current-valley, states with 
higher energies become localized on the right side of the 
source tip. So the overlap of LDOS's on the ends of the 
source tip is reduced. As a result, their product which 
is proportional to the transmission decreases. If these 
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FIG. 6: (Color online) Local density of states on the first 
(atom located on 20A° , long dashed line) and last (atom lo- 
cated on 25 A° , dashed line) atoms of the source tip is plotted 
in upper half of the graphs. Their product (dotted line) is 
compared with the transmission (solid line) in lower half of 
the graphs. Voltages are at the current-peak (0.4 V, left) 
and current-valley (0.55 V, right). The chain connected to 
the Graphene tips contains 4 atoms. Vertical lines show the 
Fermi level and the integration window. For comparison with 
LDOS products, transmission is shown 10 3 times larger. 
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FIG. 7: In the image charges method, test charge induces 
some image rings just behind the boundary surface. 



VII. APPENDIX 



The variational Method 



localized states fall in the integration window of current, 
transmission as well as current reduction occurs. 



To find the effect of image charges, we need to impose 
the Dirichlet boundary condition V — at the two left 
and right electrode planes, instead of solving Poisson's 
equation, we postulate the electrostatic Green's function 
of Eq.flTJ to be: 



VI. CONCLUSION 



In conclusion, for observation of NDR, although the 
presence of sharp features in the density of states located 
on the sharp tip apexes and their localization is required, 
the enhancing factor for localization is the charge deple- 
tion of the molecule as the bias is increased. Asymmetric 
potential profile which shows a sharp potential drop in 
the source side of the molecule, arises from the asymme- 
try in the LDOS of electrodes connected to the molecule. 
The asymmetry in LDOS's causes different amounts of 
charge flow from the molecule to the drain and source 
electrodes, respectively. The weak screening of the po- 
tential due to the depleted charge causes a larger poten- 
tial drop on the source side. However, the potential on 
the drain side varies weakly and remains almost flat. Be- 
cause of the potential drop in the source tip, states with 
higher energy become localized on the sites with higher 
potentials (right side of the source tip), and states with 
lower energy become localized on the sites with lower po- 
tentials (left side of the source tip), similar to an ionic 
bond. The charge depletion and potential drop are inten- 
sified in the source tip as the applied voltage is increased. 
This results in a more effective localization of states. Lo- 
calization causes a reduction in the overlap of the LDOS's 
on the ends of the source tip and a subsequent reduction 
in the transmission and current. 



U{rR,rj) + U(r L ,rj) 



U{?i,?j, 2 

U(n,rj) -U(f L ,rj) 



Zi ^> Zj 



Zi <C Zn 



where Fr and ?l show the positions of the atomic lay- 
ers located in the right and left contact surfaces, respec- 
tively. Although this function is not the exact solution 
of Poisson's equation, it has the correct limits for on 
the boundary surfaces, where it is equal to zero by con- 
struction. It is therefore a reasonable solution in a varia- 
tional sense, though here we are not varying any parame- 
ter to optimize the solution. In this method, we postulate 
that the image charges potential on the test charge plane 
(zi = Zj) to be as an interpolation of the left and right 
solutions in Eq. (TT5|) . The kernel used for the coulombic 
function U has been chosen to be as the OK model in 
Eq.©. 



B. Numerical Method of Images 

The straightforward way for providing an electrostatic 
Green's function which satisfies Dirichlet boundary con- 
dition, is to use image charges. Image charges can be 
put on fictitious planes just behind the plane on which 
we want the potential to be zero. Note that the choice 
of their location or charge is not unique. 

Since the potential on the boundary surfaces must be 
zero, one can find the image charges, if their location is 
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Site Number 

FIG. 8: Comparison of three methods for the calculation of 
the electrostatic Green function. The sample is a (5,5) nan- 
otube which has 4 rings (40 atoms) in the middle part. The 
test charge is set on the site number 19. The Hubbard term 
is considered to be Uh = 11.3. Numerical calculation of the 
image charges method has been done by n = 20 and zo = 2. 

fixed, by solving a system of the linear equations. For a 
test charge located on a molecular site fj , one has to solve 
the set of linear equations which are equal in number to 
the number of boundary constraints. The constraints 
leading to a linear system are as follows: 

Vin,^) = U^r^+^tiWiA) (20) 
fe=i 



V(f L ,fj) = V(f R ,f J ) = (21) 



where f is the field point and fj is the source point, 
with its images being of charge q\ and located at p|. For 
a given test charge location, the number of images ni mg 
we need depends on the number of points (constraints) 
on the boundaries, at which one wants the potential to 
be zero. 

As an example, Fig. (Q shows a nanotube and the po- 
sition of its contacts and image charges. In this model, 
all image rings are placed behind the first image plane 
marked by number 1. The first image charge planes 
which are the reflected planes from the contact surfaces, 
are located at z — — d and z = 2L — d, where d is the 
distance of the plane which includes the test charge from 
the left contact surface. The distance of image planes 
from each other is considered to be a constant value Zq. 
The number of image planes is equal to the number of 
boundary rings (n). It is supposed that the number of 
sites on an image ring is the same as the boundaries and 
nanotube rings. In this case, cylindrical symmetry of the 
images and boundaries sites is important to produce a 
smooth potential at the boundaries. 

Fig.® shows a comparison between these three meth- 
ods. A good correspondence can be observed between the 
potential of image charges method and the exact method. 
They differ by only 2 percent, while they have about 20 
percent difference with the variational method. However, 
the advantage of the variational method is its simplicity 
for application on any structure, while the position and 
values of image charges depend on the structural symme- 
tries. 
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